home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / fractals / rm2.lha / testpoint.asm < prev    next >
Encoding:
Assembly Source File  |  1994-01-01  |  6.5 KB  |  242 lines
  1. ; SCDB:
  2. ; The easiest way to build an assembly language module is to let the compiler
  3. ; do it, and hand-optimize from there. This program is a particularly dramatic
  4. ; case of the 90/10 rule: 90% of the time is spent in 10% of the code.
  5. ; In this program, it's more like 99/1. However, we pay a
  6. ; price for my laziness in starting with the compiler output, and
  7. ; that is that the branch points have cryptic names line ".10004".
  8. ; C'est la vie. (I bet I didn't spell that right. I studied German,
  9. ; not French.)
  10. ;
  11. ; The Aztec compiler actually has a pretty good optimizer for normal code in
  12. ; the regular CPU. However, the native floating point code it created was
  13. ; MISERABLE. The Aztec compiler spent excessive amounts of time moving values
  14. ; onto and off of the stack, and only used fp0. There are actually 8
  15. ; floating point registers, which means we can use most of them as
  16. ; "register double" variables, not competing with pointer or integer variables
  17. ; in the regular CPU. (I have to consider this a major flaw in the compiler.)
  18. ; Manually, therefore, I do exactly that as follows:
  19. ;
  20. ;    fp7 = x1
  21. ;    fp6 = y1
  22. ;    fp5 = x2
  23. ;    fp4 = y2
  24. ;    fp3 = x_squared
  25. ;    fp2 = y_squared
  26. ;
  27. ; The values "cx" and "cy" still come off the stack, but they should be
  28. ; cache-hits, so they ought to run pretty fast.
  29. ;
  30. ; fp0 is used for certain intermediate calculations. fp1 contains the value
  31. ; above which the loop terminates. For a standard Mandelbrot diagram, this value
  32. ; is 4.0, but it occurs to me that that value is rather arbitrary, so
  33. ; I'm putting it into fp1 so that we can play games with it. Just what happens
  34. ; if that value is changed? Beats heck out of me, but it might be interestng
  35. ; to find out... This is therefore set up to load this value from
  36. ; a double variable from the C code.
  37. ;
  38. ; To the extent that I understand it, I'm also trying to optimize this
  39. ; to take advantage of the 68882 concurrency capability. This is completely
  40. ; compatible with a 68881 - it just runs less rapidly since the '881 has no
  41. ; such concurrency capability. (On the '882, certain FMOVE.D operations can
  42. ; take place simultaneously with real operations and therefore take no
  43. ; excess time.)
  44. ;
  45. ; On a 68040, this should all run out of cache, which will speed it up still
  46. ; further. It may or may not fit into the 68030 cache. Someone out there
  47. ; who owns an '030 may wish to pursue this.
  48. ;
  49. ; To cut down on excessive FMOVE.D commands, the loop is unrolled into a double
  50. ; cycle (as per recommendations in the 68882 technical manual from Motorola).
  51. ; This costs a miniscule amount of space but saves a miniscule amount of time
  52. ; per loop, multiplied by millions of loops per picture. I consider this a
  53. ; very good trade. If I can trade a kilobyte for 5 seconds, that's a good deal,
  54. ; and the improvement is considerably better than that for some pictures.
  55. ;
  56. ; This code is in format for the Aztec assembler. Your assembler may vary
  57. ; considerably from this. The C-code given in comments below this point does
  58. ; NOT match that from the program I began with. What you're looking at is
  59. ; the C-code I ended up with after I did as much optimization as I could at
  60. ; the C-code level (including unwrapping the loop, and eliminating redundant
  61. ; calculation of "x_squared" and "y_squared").
  62. ;
  63. ; Code lines in lower case are retained from Aztec. Upper case is mine.
  64. ;
  65. ;:ts=8
  66.     far    code
  67.     far    data
  68.     machine    mc68020
  69.     mc68881
  70. ;#include <math.h>
  71. ;extern short max_orbits;
  72. ;extern double threshold;
  73. ;
  74. ;/* Test_Point() tests points in the C plane to determine
  75. ;   whether they orbit stably or escape to infinity.
  76. ;   Returns the escape velocity or -1 if the point is
  77. ;   in the M-set. */
  78. ;
  79. ;short int Test_Point(double cx, double cy)
  80. ;{
  81.     xdef    _Test_Point
  82. _Test_Point:
  83.     link    a5,#.2
  84.     movem.l    d2,-(sp)
  85.  
  86. ; SCDB: By definition, "x" and "y" begin at zero at the first loop.
  87. ; Therefore, by definition, "x_squared" and "y_squared" also start at 0.
  88. ;
  89. ;  double x1 = 0.0, y1 = 0.0;
  90. ;  double x2, y2;
  91. ;  register short int i;
  92. ;  double x_squared = 0.0, y_squared = 0.0;
  93. ;
  94. ;  i = max_orbits;
  95. ;
  96. ; SCDB: Aztec makes "int" a long, but for our purposes we'll assume it is
  97. ; a short. The only case in which this is a problem is if "/t" in the regular
  98. ; code exceeds 126, in which case the picture will probably take a year to
  99. ; calculate anyway. The following move.w therefore discards the top 16 bits.
  100.     move.w    _max_orbits+2,d2
  101.  
  102. ; SCDB: Load the threshold, nominally "4.0" into 'fp1' for future use:
  103.  
  104.     FMOVE.D    _threshold,fp1
  105.  
  106. ;
  107. ; The Aztec compiler says that the following constant is "0.0". Who am I to
  108. ; doubt it?
  109.  
  110.     FMOVE.D    #"$0000000000000000",fp7    ; x1 = 0
  111.  
  112.     FMOVE.X    fp7,fp6        ; y1 = 0
  113.     FMOVE.X    fp7,fp3        ; x_squared = 0
  114.     FMOVE.X    fp7,fp2        ; y_squared = 0
  115.  
  116.  
  117. ;  for (;;) {
  118. ;
  119. ;    y2 = x1*y1;
  120. .10003
  121.  
  122.     FMOVE.X    fp7,fp4        ; y2 = x1
  123.     FMUL.X    fp6,fp4        ;         * y1
  124.  
  125. ;    y2 = y2 + y2 + cy;
  126.  
  127.     FADD.X    fp4,fp4        ; y2 = y2 + y2
  128.     FADD.D    16(A5),fp4    ;              + cy
  129.  
  130. ;    x2 = x_squared - y_squared + cx;
  131.  
  132.     FMOVE.X    fp3,fp5        ; x2 = x_squared
  133.     FSUB.X    fp2,fp5        ;                - y_squared
  134.     FADD.D    8(A5),fp5    ;                            + cx
  135.  
  136. ;    if ((x_squared=(x2*x2)) + (y_squared=(y2*y2)) > threshold) {
  137.  
  138.     FMOVE.X    fp5,fp3        ; x_squared = x2
  139.     FMUL.X    fp5,fp3        ;                * x2
  140.  
  141.     FMOVE.X    fp4,fp2        ; y_squared = y2
  142.     FMUL.X    fp4,fp2        ;                * y2
  143.  
  144.     FMOVE.X    fp3,fp0        ; temp = x_squared
  145.     FADD.X    fp2,fp0        ;                  + y_squared
  146.  
  147.     fcmp.x    fp1,fp0
  148.     fble    .10004        ; Branch if we're NOT done
  149.  
  150. ; fall through if we ARE done
  151.  
  152. ;      return(max_orbits - i);
  153.  
  154.     move.l    _max_orbits,d0
  155.     move.w    d2,a0
  156.     sub.l    a0,d0
  157.  
  158. .3
  159.     movem.l    (sp)+,d2
  160.     unlk    a5
  161.     rts
  162. ;    }
  163.  
  164.  
  165. ;    if ((--i)<=0)
  166. ;    return(-1);
  167. .10004
  168.     sub.w    #1,d2
  169.  
  170.     bgt    .10005
  171.  
  172.     move.l    #-1,d0
  173.     bra    .3
  174.  
  175.  
  176. ;
  177. ;    y1 = x2*y2
  178. ;
  179.  
  180. .10005 
  181.  
  182.     FMOVE.X    fp5,fp6        ;y1 = x2
  183.     FMUL.X    fp4,fp6        ;        * y2
  184.  
  185.  
  186. ;    y1 = y1 + y1 + cy;
  187. ;
  188.  
  189.     FADD.X    fp6,fp6        ; y1 = y1 + y1
  190.     FADD.D    16(a5),fp6    ;              + cy
  191.  
  192. ;    x1 = x_squared - y_squared + cx;
  193.  
  194.     FMOVE.X    fp3,fp7        ; x1 = x_squared
  195.     FSUB.X    fp2,fp7        ;                - y_squared
  196.     FADD.D    8(a5),fp7    ;                            + cx
  197.  
  198. ;    if ((x_squared=(x1*x1)) + (y_squared=(y1*y1)) > threshold) {
  199.  
  200.     FMOVE.X    fp7,fp3        ; x_squared = x1
  201.     FMUL.X    fp7,fp3        ;                * x1
  202.  
  203.     FMOVE.X    fp6,fp2        ; y_squared = y1
  204.     FMUL.X    fp6,fp2        ;                * y1
  205.  
  206.     FMOVE.X    fp3,fp0        ; temp = x_squared
  207.     FADD.X    fp2,fp0        ;                  + y_squared
  208.  
  209.     FCMP.X    fp1,fp0
  210.     fble    .10006
  211.  
  212. ;      return(max_orbits - i);
  213.  
  214.     move.l    _max_orbits,d0
  215.     move.w    d2,a0
  216.     sub.l    a0,d0
  217.     bra    .3
  218. ;    }
  219.  
  220. .10006
  221. ;    if ((--i)<=0)
  222. ;    return(-1);
  223.  
  224.     sub.w    #1,d2
  225.  
  226.     bgt    .10007
  227.     move.l    #-1,d0
  228.     bra    .3
  229. ;  }
  230. .10007
  231. .10001
  232.     bra    .10003
  233. ;}
  234. .2    equ    -48
  235. ;
  236.     xref    .begin
  237.     dseg
  238.     xref    _max_orbits
  239.     xref    _threshold
  240.     end
  241.  
  242.